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Abstract 


High-order  essentially  non-oscillatory  (ENO)  finite-difference  schemes  are  applied  to  the 
two-  and  three-dimensional  compressible  Euler  and  Navier-Stokes  equations.  Practical  issues, 
such  as  vectorization,  efficiency  of  coding,  cost  comparison  with  other  numerical  methods 
and  accuracy  degeneracy  effects,  are  discussed.  Numerical  examples  are  provided  which  are 
representative  of  computational  problems  of  current  interest  in  transition  and  turbulence 
physics.  These  require  both  non-oscillatory  shock  capturing  and  high  resolution  for  detailed 
structures  in  the  smooth  regions  and  demonstrate  the  advantage  of  ENO  schemes. 
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1  Introduction 


In  the  computation  of  inviscid,  compressible  flow,  the  presence  of  infinitesimally  thin  shocks 
readily  leads  to  non-linear  instability  for  traditional,  unadulterated,  linearly  stable  high- 
order  methods.  Moreover,  regions  of  strong  gradients  which  have  finite  thickness  but  are 
too  thin  for  the  grid  to  resolve  may  also  produce  non-linear  instability.  This  is  the  case, 
for  example,  for  high  Reynolds  number  Navier-Stokes  computations  in  which  the  shock 
thickness  is  much  smaller  than  the  grid  spacing.  The  standard  “cures”  are  either  to  add 
explicit  artificial  viscosity  to  the  numerical  method  or  to  employ  an  upwind-bieused  scheme 
(which  conteiins  implicit  artificial  viscosity).  Such  approaches  usually  have  the  undesirable 
side  effect  of  loss  of  resolution,  particularly  for  the  small-scale  structures  in  the  smooth  part 
of  the  solution.  The  small-scale  features  are  typically  strongly  and  erroneously  dampedhy  the 
jurtificial  viscosity.  This  problem  even  afflicts  the  formally  high-order  TVD  (total- variation- 
diminishing)  schemes,  since  they  must  degenerate  to  first  order  accuracy  at  smooth  critical 
points  [13].  Certainly,  the  most  difficult  feature  to  capture  is  the  passage  of  small-scale 
features  through  shook  waves. 

In  many  applications,  such  as  typical  steady-state  aerodynamic  CFD,  the  side  effects  of 
the  artificial  viscosity  are  not  particularly  worrisome  since  the  main  target  of  the  computa¬ 
tion  is  the  large-scale  flow  structure  and  the  details  of  the  flow  near  the  shock  are  not  too 
significant.  This  is  decidedly  not  the  case,  however,  for  numerical  simulations  of  transition 
and  turbulence.  Here  the  interesting  physical  phenomena  occur  on  scales  much  smaller  than 
those  of  the  mean  flow.  As  noted  by  Hussaini  and  Zang  [9],  for  incompressible  flow  spectral 
methods  have  been  preferred  for  these  applications  since  they  have  the  best  fidelity  for  the 
small-scale  flow  features.  However,  due  to  their  extreme  sensitivity  to  non-linear  instability 
spectral  methods  have  yet  to  be  used  for  serious  investigations  of  transition  and  turbulence 
in  compressible  flow  with  regions  of  strong  gradients.  (They  can  be  used,  of  course,  for 
shock-free  flows  [4]  and  for  low  Reynolds  number  viscous  flows  [5]  in  which  thick  shocks  are 
actually  resolved  rather  than  captured.) 

Essentially  non-oscillatory  (ENO)  schemes,  first  introduced  by  Harten  and  Osher  [6]  and 
Harten,  Engquist,  Osher  and  Chakravarthy  [7],  can  achieve  uniformly  high-order  accuracy 
with  sharp,  essentially  non-oscillatory  shock  transitions.  The  key  idea  is  an  adaptive  sten¬ 
cil  interpolation  (based  on  difference  tables)  which  automatically  interpolates  in  a  locally 
smoothest  region.  This  strategy  provides  a  strong  inhibition  towards  differencing  across 
discontinuities.  In  [22,  23],  Shu  and  Osher  introduced  an  efficient  implementation  of  ENO 
schemes,  using  the  same  adaptive  stencil  idea  but  working  directly  on  fluxes  and  a  special 
class  of  TVD  high-order  Runge-Kutta  type  time  discretizations.  It  bypasses  the  reconstruc¬ 
tion  and  Lax-Wendroff  procedures  in  the  original  ENO  schemes.  For  multi- dimensions,  this 
simplification  is  significant,  because  the  reconstruction  in  multi-dimensions  becomes  quite 
complicated  [8].  Numerical  examples  in  [22]  and  [23],  especially  the  examples  of  shock  in¬ 
teraction  with  entropy  and  vorticity  waves,  for  which  a  good  resolution  for  the  detailed 
structures  in  the  smooth  region  is  as  important  as  a  sharp,  non-oscillatory  shock  transition. 
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indicate  a  good  potenti;d  for  ENO  schemes  in  computing  compressible  Euler  and  Navier- 
Stokes  equations. 

In  this  paper,  we  discuss  the  coding  of  ENO  schemes  in  [22,  23]  to  two  dimensional  general 
geometry  (via  transformation)  arid  to  three-dimensional  Euler  and  Navier-Stokes  equations 
of  compressible  gas  dynamics  on  Cray  X-MP  and  Y-MP.  We  address  the  practical  issues  such 
as  vectorization,  efficiency  of  evaluating  Newton  interpolations,  cost  comparison  with  other 
numerical  methods,  and  accuracy  dego'ieracy  effects.  We  then  present  numerical  examples 
which  all  require  ncn-oscillatory  shock  capturing  and  high  resolution  for  rich  structures 
in  the  smooth  regions,  including  two-dimensional  shear  flows,  two  dimensional  and  three 
dimensional  homogeneous  turbulence,  aid  two-dimensional  shock  interaction  with  entropy 
and  vorticity  waves. 


2  The  Navier-Stokes  and  Euler  equations 


For  completeness,  we  document  here  the  three-dimensional,  compressible  Navier-Stokes  equa¬ 
tions  as  well  as  the  various  transformations  that  are  required  for  the  ENO  method  in  curvi¬ 
linear  coordinates.  In  terms  of  the  density  p,  the  velocity  u  =  (u,?;,u;)*,  the  pressure  p,  and 
the  internal  energy  £?,  the  Navier-Stokes  equations  read 

qt  +  ^q)*  +  g(q)v  +  h(q)*  =  r(q)«  +  s(q)v  +  t(q),  (1) 

where 

q  =  (p,pu,pu,;ic,xi;%  f(q)  =uq  +  p(0,l,0,0,u)‘ 
g(q)  =uq  +  p(0,0,l,0,u)‘,  h(q)  =ii;q  +  p(0,0,0,l,u;)‘  (2) 

and 


r(q)  =  ;i(0,rn,T2i,T3i,c.-i)*, 

s(q)  =  /i(0,Ti2,T22,T32,cr2)‘, 

t(q)  =  /i(0,Ti3,T23,T33,O-3)‘,  (3) 

with  the  components  of  the  viscous  stress  tensor  given  by 
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Also,  fi  is  the  viscosity,  7  is  the  ratio  of  specific  heats,  Pr  is  the  Prandtl  number,  and 
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For  implementing  ENO  schemes  with  characteristic  decompositions,  we  need  the  expressions 

di  dz.  5h 

of  the  eigenvalues  and  the  right  and  left  eigenvectors  of  the  Jacobians  — ,  — ,  — .  The 

oq  aq  aq 

dt 

eigenvalues  for  are  u  —  c,u,u,u,u  +  c.  Its  right  eigenvectors  are  the  columns  of 
dq 
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and  the  left  eigenvectors  are  the  rows  of 
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The  corresponding  expressions 


dg  d\i 

for  and  —  are  apparent, 
aq  aq 
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The  Euler  equations  can  be  obtained  from  the  Navier-Stokes  equations  by  setting  fi  =  0. 
The  equations  for  two-dimensional  problems  are  equally  obvious. 

In  the  two-dimensional  case,  the  transformation 

®  =  y=^yiLv)>  (10) 

enables  us  to  treat  non-uniform  grids  or  mappings  into  non-rectangular  domains.  The  Navier- 
Stokes  equations  become 

4  +  +  g„  =  +  s„,  (11) 

where 

q  =  J-% 

f  =  J-^m  +  piO,U^y,U)% 
g  =  J-^[Vq  +  p{0,rj,,ny,V)% 
r  =  J-'(^»r+4s), 
s  =  J-'^{rjJ+r)ys), 
r  =  /i(0,ru,T2i,cri)‘, 

S  =  /i(0,Ti2,T22,£r2)‘,  (12) 

with 

Tu  =  ^{^xUi  +  VxU„)-^{^yV^+VyVr,), 

T2i  =  ri2  =  +  T)yU„  -f 

4  2 

•^22  =  +  Vy-^v)  -  +  VxUr,), 

(Ti  =  urn  +  vTi2  +  J^ZIifPr  +  Vx{c^)r,], 

<T2  =  UT21  4-  VT33  +  +  Vy{c^)r,]>  (13) 

and 


J  —  CxVy  ^yVx} 

U  =  ^*U  +  CyV, 

V  =  rjxU  +  TjyV.  (14) 


The  eigenvalues  for 
the  columns  of 
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are  Z7  -  c  U,U,U  +  c  +  ^y;  the  right  eigenvectors  are 
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the  left  eigenvectors  are  the  rows  of 
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We  agadn  forego  the  explicit  prescription  for 

dq 


\ 


-2bi  '  ’ 

(62  -  -(fciu  -  -(M  -  ^)  bi  J 

where  bi  and  ^2  are  defined  in  (2.8);  and 


3  Implementing  the  ENO  schemes 


(16) 


(17) 


This  section  should  be  read  in  conjunction  with  [22, 23]  for  notation,  terminology  and  other 
details  of  ENO  schemes  based  on  fluxes  and  for  TVD  time-discretizations  of  Runge-Kutta- 
type.  Here  we  only  summarize  several  key  steps  of  the  algorithms  and  address  practical  issues 
such  as  vectorizatioa,  efficiency,  cost  comparison  and  the  reduction  of  round-off  errors. 

The  ENO  procedure  is  applied  only  to  the  convection  part,  i.e.,  the  left-hand-side,  of 
Eq  (1).  The  diffusion  righthand  side  of  Eq.  (1)  is  approximated  by  the  standard  centered 
differences.  It  is  also  possible  to  use  ENO- type  adaptive  stencil  interpolation  to  approximate 
the  diffusion  terms,  but  we  have  not  observed  any  significant  differences  in  our  numerical 
tests  (typically  with  small  physical  viscosity  u).  Besides  simplicity,  centered  approximations 
also  seem  more  natural  for  diffusion  terms. 


We  now  summarize  the  key  steps  of  the  algorithm; 

(1)  The  time  marching  is  implemented  by  a  class  of  TVD  Runge-Kutta  type  methods  [22]. 
For  example,  the  third  order  case  is 

^(1)  ^  +  AtL(q^°^), 

e-  =  |q‘”'  +  iqW  +  iA(L(g''>), 
qW  =  i^°'  +  fq'“>  +  5AtL(qW), 

5"+'  =  (18) 


where  L  is  the  numerical  spatial  operator  approximating  the  spatial  derivatives  in  Eq.  (2). 
This  class  of  Runge-Kutta  methods  is  labeled  TVD  because  it  has  been  shown  [22]  that  it 
does  not  increase  the  total  variation  of  the  spatial  part  under  a  suitable  CFL  restriction. 
Also  notice  that  for  the  third  order  case,  Eq.  (18),  only  three  storage  levels  (two  for  q,  one 
for  L)  are  needed,  since  can  overwrite  q^^^  and  q^®^  can  again  overwrite  q^^^ 

(2)  We  thus  only  need  to  consider  the  spatial  operator 

L(q)  =  -f(q)i  -  g(q)„  -  h(q);,  +  r(q)*  +  s(q)„  +  f(q)*.  (19) 


The  last  three  terms  are  approximated  by  standard  second-order  or  fourth-order  centered 
differences.  We  use  an  ENO  scheme  based  on  fluxes;  hence  the  first  three  terms  can  be 
approximated  dimension-by-dimension:  for  example,  when  approximating  — f(q)*,  y  and  z 
are  fixed.  The  core  of  the  algorithm  is  then  a  one-dimensional  ENO  approximation  to,  e.g. 
-f(q)x. 

(3)  Since  f (q)  is  a  vector,  we  can  approximate  — f(q)*  either  component  by  component, 
or  in  (local)  characteristic  directions.  In  the  former  case,  to  obtain  non-linear  stability  by 


upwinding,  we  write 

?(q)  =  f^(q)  +  /“(q))  (20) 

with 

<^(q)  =  2(^(^)  (21) 

d{ 

where  a  =  max(|u|  -f-  c)  is  the  largest  eigenvalue  in  absolute  value  of  the  Jacobian  ^ 

d{^ 

along  the  relevant  s-line.  The  decomposition  in  Eq.  (20)  guarantees  that  has  posi¬ 


tive/negative  eigenvalues;  hence  upwinding  (to  be  discussed  later  in  this  section)  is  the  same 
for  all  components.  For  characteristic  decompositions,  we  take  Aj+1/3  to  be  some  average 
Jacobian  at  Sj+1/2,  e.g.  the  arithmetic  mean 


^^+5  -  dq 


q=5(q,+C+i) 


or 


- 


dq 


_  -.Roe  ’ 

q=q.+i 


(22) 


(23) 


.(Roe) 


where  q .  i  ^  is  the  Roe  average  of  and  [17].  We  then  use  the  left  eigenvectors 

Rjli/2  ill  Eq.  (8)  or  Eq.  (16)  of  Aj+1/2  to  project  all  relevant  quantities  (differences  around 
Xj+1/2)  to  the  local  characteristic  fields.  A  scalar  ENO  algorithm  can  then  be  applied,  and 
the  result  projected  back  to  component  space  by  J2j+i/2  in  Eq.  (7)  or  Eq.  (16). 
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(4)  We  finally  describe  the  implementation  of  the  scalar  ENO  approximation  of  — /(?)*. 
It  is  written  as  a  conserved  flux  difference 

-f{<l)x\x=xj  ^  ?+5  ~  ^ 3-0' 

where  the  numerical  flux  fj+1/2  approximates  h{xj+y2)  to  a  high  order  with  h{x)  defined  by 

It  is  pointed  out  in  [23]  that  we  do  not  need  to  construct  1  (a)  explicitly:  we  simply  use  the 
difference  tables  of  f(u{x)).  If  the  (undivided)  differences  cf  f(u{x))  are  defined  by 


/b'.O]  =  f(uj), 

f[3,k]  =  f\j  +  l,k-l]-  f\j,k-l], 


where  r  is  the  spatial  order,  then 


fj+i  =  X)  c(i-i,m)/[i,m], 


with  i  being  the  left-most  point  in  the  stencil  used  to  approximate  fj+i/2)  c(s,m)  being 
defined  by 

1  *+m  •+m 

=  n  (-?)•  (28) 

*“*  “P  —  3 
P  7^  I 

The  small,  constant  matrix  c  is  computed  once  and  stored. 

The  adaptive  stencil  determined  by  the  choice «,  +  the  left-most  point  in  the  stencil  used 
to  approximate  fj+i/2-  We  start  with  i  =  j  or  i  =  j-\  l  according  to  the  (local)  wind 
direction  (upwinding),  and  then  apply  the  following 

if  (absCf [i,k]) .gt.absCf Ci-l,k]))  i=i-l  (29) 


for  A:  =  1,  -  •  •  ,r. 

Remark  3.1  The  code  is  written  in  such  a  way  that  all  the  major  loops  are  vectorized  by 
default  of  Cray  Fortran.  To  vectorize  Eqs.  (27)-(29)  we  can  either  repeat  (29)  r  times  (for 
fixed  r  only)  or  introduce  a  temporary  one  dimensional  storage  for  i  to  put  the  short  loop  (29) 
outside  the  long  loop  (27)-(29)  for  j.  To  vectorize  the  characteristic  decompositions  we  have 
to  introduce  one  dimensional  temporary  storage  for  the  local  projection  on  characteristic 
fields  at  each  xj.  Since  we  only  vectorize  the  innermost  one-dimensional  loop,  we  need  just 
17  three  dimensional  units  (10  for  two  components  of  q,  5  for  L,  and  2  work  units)  for  three 
dimensional  problems  using  third  order  schemes. 


Remark  3.2  For  our  current  implementation,  the  componentwise  ENO  scheme  takes 
around  4.5  times  as  much  CPU  time  as  a  centered  finite-difference  scheme  with  the  same 
order  of  accuracy.  A  factor  of  two  is  due  to  the  flux  splitting,  Eq.  (20).  Instead  of  just 
computing  f  (q)i  we  are  computing  t  '  (q)*  +  f  (q)i;  hence  the  work  is  doubled.  This  is 
the  price  to  pay  for  implementing  upwinding  techniques  to  achieve  stability.  Another  factor 
of  two  is  due  to  the  adaptive  stencil  procedure,  Eq.  (29):  when  these  “if”  statements  are 
removed,  the  code  runs  twice  as  fast.  It  seems  odd  that  these  “if”  statements  account  for  so 
much  CPU  time  since  they  are  all  vectorized.  The  main  reason  is  that  since  the  pattern  is  not 
uniform  from  point  to  point,  gathering  and  scattering  are  activated  by  Cray  Fortran.  These 
procedures  are  very  slow  on  the  Cray.  A  similar  slow-down  also  exists  for  TVD  schemes. 
ENO  schemes  using  characteristic  decompositions  take  more  CPU  time:  ENO-LF  and  ENO- 
Roe  with  entropy  fix  (see  [23]  for  definitions)  take  about  3  and  1.7  times,  respectively,  as 
much  CPU  time  as  componentwise  ENO  .schemes.  Notice  that  ENO-Roe  is  faster  than 
ENO-LF  because  it  does  not  use  a  flux  splitting.  See  [23]  for  more  detauls. 

Remark  3.3  We  use  undivided  differences,  Eq.  (26),  and  prestored  local  matrix  c,  Eq. 
(28),  to  reduce  cost  and  to  reduce  the  effect  of  round-off  errors. 


4  Transition  in  a  Free  Shear  Layer 


The  numerical  examples  were  chosen  to  illustrate  the  ENO  method  for  problems  in  transition 
and  turbulence.  We  consider  flows  with  gradients  which  are  readily  resolved  by  the  grid  - 
shocks  are  either  absent  altogether  or  else  sufficiently  broad  (due  to  low  Reynolds  number) 
that  the  usual  spectral  method  gives  stable  results.  Comparisons  of  the  intrinsic  resolution 
can  therefore  be  made  between  the  spectral  and  ENO  methods..  We  then  take  strong  shock 
cases  to  illustrate  the  advantage  of  non-oscillatory  high-order  methods. 

Unless  otherwise  indicated,  third-order  ENO  with  the  third-order  Runge-Kutta  time  dis¬ 
cretization  (18)  and  fourth-order  centered  differences  for  t.ie  viscous  terms  are  used.  Notice 
that  the  third-order  ENO  [23]  is  actually  fourth  order  in  smooth,  monotone  regions;  hence 
for  problems  with  isolated  critical  points  it  is  fourth-order  in  Lx  norm  sense.  We  use,  as 
in  [23],  the  notations  ENO-LF  (Lax-Friedrichs),  ENO-Roe  and  ENO-Com  (componentwise). 

There  has  been  considerable  recent  interest  in  the  physics  of  the  compressible  free  shear 
layer  and  numerical  simulations  have  furnished  several  interesting  results.  The  numerical 
methods  employed  have  typically  been  second-order  TVD  [24,  12],  fourth-order  MacCor- 
niack  [25,  16],  or  fourth-  or  even  sixth-order  compact  [10,  19,  2].  Atkins  [1]  and  Sandham 
and  Yee  [20]  have  made  detailed  studies  of  the  performance  of  TVD  schemes  on  this  prob¬ 
lem.  The  latter  also  made  comparisions  with  sccond-ordcr  MacCormack  results.  Carpenter, 
et  al.  [2]  have  compared  third-order  upwind,  fourth-order  MacCormack  and  fourth-order 
compact  methods. 

For  the  particular  2-D  examples  studied  in  the  present  paper,  the  mean  flow  is  given  by 
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the  hyperbolic  tangent  profile  Uq  =  tanh(y),  vo  =  0,  and  po  =  where  Moo  is  the  ratio 
in  the  limit  y  —*  ±oo,  of  the  free  stream  velocity  to  the  sound  speed,  and  periodic  boundary 
conditions  are  enforced  in  the  streamwise  (x)  direction.  The  velocity  is  non-dimensionalized 

duo 

by  the  freestream  velocity,  Uoo,  lengths  by  half  the  vorticity  thickness  6^  =  2uo/|-j— |,  the 

oy 

density  by  the  freestream  density,  poo,  the  temperature  by  the  freestream  temperature.  Too, 
and  the  pressure  by  Poo^iL-  The  Reynolds  number  Re  =  uqS^Poo! p^oo-  The  viscosity  p.  is 
prescribed  through  Sutherland’s  law  with  a  reference  temperature  of  520‘’X  and  the  Prandtl 
number  is  taken  to  be  0.7.  This  example  is  for  the  temporally  evolving  free  shear  layer.  A 
forcing  term  is  added  to  the  Navier- Stokes  equations  in  order  to  make  the  assumed  mean 
flow  a  steady  solution.  The  computational  domain  is  (0,27r/a)  x  (—00,00),  where  ot  is  a 
specified  wavelength. 

For  this  problem  we  present  comparisons  of  ENO  with  both  explicit  and  compact  centered 
difference  schemes  and  with  spectral  methods.  The  explicit  central  difference  methods  use 
3  points  for  a  second-order  approximation  and  5  points  for  fourth-order.  The  compact 
difference  scheme  uses  a  Fade  approximation  with  5  explicit  points  and  3  implicit  points 
(see  [10]).  It  is  formally  sixth-order  accurate  at  all  interior  points,  and  at  points  adjacent  to 
the  boundary,  but  reduces  to  fourth-order  accuracy  at  the  boundary  itself.  (The  sixth-order 
compact  scheme  used  by  Lele  [10]  reduces  to  fourth-order  accuracy  at  points  adjacent  to 
the  boundary  and  third-order  accuracy  at  the  boundary  itself.)  For  the  spectral  calculation 
a  Fourier  expansion  is  applied  in  x  and  the  Cain  transformation  [3] 

y  =  —L  cot(77)  (30) 

is  used  in  the  y  direction;  L  is  a  stretching  parameter  which  is  taken  between  4  and  10. 
This  permits  the  use  of  cosine  (for  p,  u,  and  e)  and  sine  (for  u)  expansions  of  the  dependent 
variables  in  the  r]  direction. 


4.1  Linear  Instability 

For  the  smooth  problem  we  consider  first  the  evolution  of  a  small  perturbation,  with  stream- 
wise  wavenumber  a  =  0.4,  from  the  mean  flow  at  Moo  =  0.5  with  a  Reynolds  number  of 
100,  and  the  usual  7  =  14.  The  shape  of  the  perturbation  is  given  by  the  fastest-growing 
eigenfunction  of  the  linearized  Navier-Stokes  equations  at  the  specified  wavelength.  The 
amplitude  of  the  perturbation  was  chosen  so  that  its  transverse  velocity  at  y  =  0  was  0.1% 
of  the  freestream  velocity.  The  growth  rate  for  this  particular  case  was  0.127464941.  (This 
eigenvalue  problem  was  solved  by  a  spectral  linear  stability  code  [11].)  For  such  small  am¬ 
plitudes  the  non-linear  code  should  produce  linear  growth  of  the  perturbations  for  small 
times.  Table  1  shows  the  growth  rates  produced  by  central  difference  methods  and  compact 
methods  after  10  time-steps  with  a  time-step  of  0.01.  (The  measured  growth  rate  was  taken 
to  be  cr  =  [log  £^(<)  -  log  jE;(0)]/(2t)  where  E{t)  =  |lu-uo|l|,2  +  11^11^3.)  In  these  examples  a 
highly-resolved  discretization  (with  128  points)  was  employed  in  y  and  the  specified  method 
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N.  =  16 


32 


method 


2nd-order  central 
4th-order  central 
4th-order  compact 
6th-order  compact 


N^  =  8 


-4.26 

-4.92 

2.15 

1.12 


-1.24 

-1.35 

5.00 

-2 


Table  1:  Linear  Growth  Rate  Errors  for  Central-Difference  and  Compact  Schemes  at  t  =  0.1 


method 

00 

11 

Nx  =  32 

Nx  =  64 

2nd-order  ENO 
2nd-order  ENO-2 

.1.14(-3) 

-8.28(-4) 

-4.90(-4) 

■5.58(-4) 

3rd-order  ENO 
3rd-order  ENO-2 

B 

-9.85(-5) 

-2.51(-6) 

4th- order  ENO 
4th-order  ENO-2 

-4.96(-3) 

-4.21(-3) 

-1.88(-4) 

-1.28(-5) 

-1.54(-5) 

-1.39(-5) 

-9.88(-6) 

-1.78(-5) 

Table  2;  Linear  Growth  Rate  Errors  for  ENO  Schemes  at  t  =  0.1 

was  applied  in  x  with  the  streamwise  resolution  as  noted  in  the  table.  The  table  thus  pro¬ 
vides  the  accuracy  achieved  as  a  function  of  the  method  and  the  number  of  grid-points  per 
wavelength.  (Even  for  =  4,  the  error  from  a  spectral  discretization  in  x  is  already  smaller 
than  10"^.) 

Similar  results  for  ENO  methods  of  second-,  third-,  and  fourth-order  are  given  in  Table 
2.  In  these  cases  the  ENO  method  was  also  applied  in  y,  again  with  128  points  used  in 
this  direction.  Again,  the  intent  was  to  isolate  the  discretization  errors  in  x.  Except  for 
the  Nx  =  64  cases,  the  results  are  as  one  would  expect:  the  accuracy  increases  with  the 
number  of  grid-points  and  with  the  order  of  the  scheme.  (The  unexpected  deterioration  of 
the  convergence  rate  for  the  finest  grid  is  addressed  below.)  A  comparison  of  the  fourth- 
order  central-difference  results  from  Table  1  with  those  of  the  3rd-order  ENO  from  Table  2 
indicates  that  the  ENO  results  are  slightly  less  accurate.  This  is  to  be  expected  since  one 
anticipates  that  in  most  of  the  flow  this  ENO  stencil  reduces  to  the  fourth-order  central  one, 
and  where  it  is  switched  to  one-sided  it  will  lose  one  order  in  accuracy. 

Notice  that  the  improvement  in  accuracy  obtained  from  going  from  32  to  64  streamwise 
grid-points  is  less  than  expected.  For  the  central-difference  schemes  this  occurs  on  the  10~® 
level,  whereas  H  occurs  an  order  of  magnitude  or  more  earlier  for  the  ENO  schemes.  (This  is 
even  more  apparent  in  the  ENO  results  at  later  times,  as  evidenced  by  the  data  in  Table  .3.) 
This  is  due  to  time-differencing  and  linearization  errors  for  the  central  difference  methods, 
whereas  it  is  caused  by  the  loss  of  accuracy  when  the  stencil  switches  for  the  ENO  method. 
The  results  labelled  by  “ENO-2”  in  the  tables  are  computed  by  using  a  factor  of  2  to  multiply 
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method 

II 

oo 

K  =  16 

11 

CO 

to 

AT.  =  64 

2nd-order  ENO 

-4.26(-2) 

-6.37(-3) 

-1.26(-3) 

2nd-order  ENO-2 

-3.86(-2) 

-1.43(-3) 

-8.33(-4) 

-5.72(-4) 

3rd- order  centered 

-4.22(-2) 

-6.21(-3) 

-1.23(-3) 

-5.89(-4) 

3rd-order  ENO 

-1.59(-2) 

-L16(-3) 

-2.41(-4) 

-1.89(-4) 

3rd-order  ENO-2 

-7.31(-3) 

1.55(-5) 

-2.52(-5) 

-3.49(-5) 

4th-order  centered 

3.55(-5) 

-3.14(-5) 

-3.56(-5) 

4th- order  ENO 

-4.95(-3) 

-2.53(-4) 

-9.35(-5) 

4th-order  ENO-2 

-4.19(-3) 

-3.12(-5) 

5th-order  centered 

-4.92(-3) 

-2.14(-4) 

-4.63(-5) 

Table  3:  Linear  Growth  Rate  Errors  at  t  =  1.0 

either  the  first  or  the  second  abs  term  in  Eq.  (29),  depending  upon  whether  i  is  greater 
than  the  left-most  point  in  the  centered  stencil  or  not.  The  effect  is  to  bias  the  scheme 
towards  a  centered  stencil  in  smooth  regions.  This  modification  of  ENO  is  discussed  in 
detail  in  [21],  accompanied  by  numerical  tests  on  smooth  and  shocked  cases,  in  response  to  a 
recent  discovery  by  Rogerson  and  Meiberg  [18]  about  some  accuracy  degeneracy  phenomena 
of  ENO  schemes.  From  the  table  we  can  see  that  “ENO-2”  is  in  most  cases  comparable  in 
accuracy  with  the  corresponding  centered  schemes,  while  ENO  is  usually  one  order  lower,  as 
expected  from  the  (unnecessary)  switching  of  stencils  in  smooth  regions. 


4.2  Mach  0.5  evolution 

Next,  we  present  a  comparison  of  these  methods  for  a  fully  nonlinear  problem.  The  previous 
results  were  just  basic  calibration  tests  (for  all  the  methods).  The  real  purpose  of  numerical 
simulation  codes  is  to  explore  nonlinear  fluid  dynamics.  The  next  example,  therefore,  is 
a  simulation  of  vortex  roll-up  and  pairing  at  M^o  =  0.5.  The  initial  conditions  consist  of 
the  mean  flow  plus  two  linear  eigenfunctions:  the  fundamental  with  wavenumber  af  =  0.4 
and  amplitude  e/  =  0.01  and  its  subharmonic  with  wavenumber  a,  =  0.2  and  amplitude 
e,  =  0.0001.  (In  this  case  the  computational  domain  in  x  is  doubled  from  that  of  the 
previous  example  in  order  to  accommodate  the  subharmonic.)  The  initial  phases  (judged  by 
the  location  of  the  maximum  of  the  normal  velocity  perturbation  at  y  =  0)  were  exactly  out 
of  phase,  a  choice  which  ensures  that  vortex  pairing  will  occur. 

Figure  1  presents  the  evolution  of  the  vorticity  thickness  for  third-order  ENO  on  grids  of 
size  32^,  64^,  and  128^  and  compares  these  results  with  those  of  a  128^  spectral  calculation. 
(An  analysis  of  the  spectral  coefficients  of  the  latter  coefficients,  along  the  lines  discussed 
in  [26],  indicates  that  the  spectral  result  is  accurate  to  better  than  4  significant  digits  until 
about  t  —  125,  but  that  thereafter  its  accuracy  deteriorates  rapidly  as  the  vortex  roll-up 
produces  scales,  particularly  in  the  streamwise  direction,  that  are  too  small  for  the  grid.) 
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The  ENO  result  is  clearly  converging  to  the  proper  answer.  A  similar  comparison  is  presented 
for  the  sixth-order  compact  scheme  in  Figure  2.  The  convergence  here  is  more  impressive 
than  for  the  ENO  scheme,  but  that  is  to  be  expected  for  this  smooth  flow.  Curiously,  the 
spectral  coefficients  for  the  compact  scheme  suggests  that  it  is  the  transverse  resolution  which 
is  most  stressed  by  the  roll-up. 


Figure  3  summarizes,  for  the  third-order  ENO  method,  the  evolution  of  the  lowest  4 
Fourier  harmonics  as  represented  by  the  quantity 


where 


Ek=  r{\nk\^  +  \vk\^)W{y)dy 

J^OO 


Qkiy,  0  =  ^  ^  9(a;.  y> 


is  the  streamwise  Fourier  coefficient  of  the  variable  q  and 

Wiy) 


_  f  1  |y|  <  yc 

\  |y|  >  ye 


(31) 

(32) 


(33) 


is  a  weight  function  used  to  confine  the  region  of  integration  to  a  finite  size.  (We  used 
Vc  =  50-)  Once  again,  the  numerical  results  are  indicative  of  convergence.  On  a  32*  grid 
the  ENO  results  are  perceptibly  different  from  the  highly  resolved  results  even  for  the  k  =  1 
mode.  On  a  64*  grid  the  worst  relative  results  occur  for  k  =  Z.  This  mode  happens  to  be 
the  most  sensitive  of  the  4  to  nonlinear  interactions.  At  the  start  of  the  calculation  only  the 
k  =  1  and  k  =  2  modes  had  non-zero  amplitudes.  The  k  =  2  mode  is  initially  forced  by 
the  self-interaction  of  the  A:  =  2  mode,  which  is  the  dominant  mode  for  the  first  part  of  the 
calculation.  The  A:  =  3  mode  is  initially  forced  by  the  interaction  between  the  A:  =  1  and 
A:  =  2  modes  and  it  is  here  that  the  largest  errors  occur.  One  heartening  result  is  that  the 
A:  =  1  mode  -  the  subharmonic  -  is  tracked  reasonably  well.  Atkins  [1]  observed  that  there 
could  be  appreciable  spurious  generation  of  this  mode  by  a  second-order  TVD  method. 


Similar  data  are  provided  in  Figure  4  for  the  compact  scheme.  The  results  for  these 
low-order  modes  are  already  graphically  indistinguishable  from  the  128*  spectral  results  on 
a  64*  grid  for  the  compact  scheme.  This,  too,  is  understandable  since  the  compact  scheme 
was  shown  in  Table  1  to  have  an  accuracy  of  better  than  1  part  in  10^  with  8  points  per 
wavelength,  and  even  the  mode  A:  =  4  has  8  points  per  wave. 

We  close  the  Mach  0.5  results  v;ith  a  plot,  in  Figure  5,  of  the  pressure  contours  at 
t  =  150  for  the  ENO  method  on  various  grids  and  for  the  high- resolution  spectral  method. 
The  similarity  of  all  the  results  is  apparent. 


4.3  Mach  0.9  evolution 

The  rationale  for  the  ENO  method  rests  primarily  on  its  behavior  in  the  presence  of  shocks. 
Indeed,  typical  central-difference  and  spectral  methods  have  substantial  difficulty  for  this 
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compressible  free  layer  problem  at  freestream  Mach  numbers  above  0.70.  For  this  example 
we  choose  Moo  =  0.9  and  initial  conditions  consisting  of  the  mean  flow  plus  two  linear 
eigenfunctions:  the  fundamental  with  wavenumber  Of  =  0.3  and  amplitude  Cf  —  0.01  and 
its  subharmonic  with  wavenumber  a,  =  0.15  and  amplitude  e,  =  0.001.  Furthermore,  the 
stretching  parameter  for  the  ENO  method  is  here  chosen  to  be  L  =  10  to  provide  better 
resolution  near  the  shock  waves  which  eventually  develop. 

The  evolution  of  the  pressure  field  for  this  case  is  depicted  in  Figure  6.  These  plots 
are  taken  from  a  computation  based  on  the  sixth-order  compact  scheme.  (A  128*  grid  was 
used  from  t  =  0  to  <  =  75,  a  256  grid  from  t  =  75  to  t  =  100,  a  512  x  192  grid  from 
t  =  100  to  i  =  106.25,  a  768  x  192  grid  from  t  —  106.25  to  <  =  112.5,  and  a  1024  x  192 
grid  from  t  =  112.5  to  t  =  137.5.  Spectral  interpolation  was  employed  for  the  requisite  grid 
refinements.)  By  i  =  100  vortex  pairing  has  already  occurred.  The  vortex  is  centered  in  the 
plotting  frame  and  stagnation  points  are  located  on  the  vertical  mid-plane  at  the  streamwise 
edges  of  the  plot.  As  discussed,  for  example,  by  Lele  [10],  the  flow  expands  away  from  the 
stagnation  points  as  it  goes  around  the  vortex  and  compresses  as  it  returns  to  the  stagnation 
point.  At  sufficiently  high  Mach  numbers  and  for  sufficiently  strong  vortices  the  compression 
occurs  via  a  shock  wave.  Shortly  after  t  =  100  in  this  case  a  pair  of  shocks  develop  -  these 
are  the  so-called  “eddy  shocklets”  [24,  10]  -  and  they  grow  steadily  stronger  as  the  flow 
evolves.  These  shocks  are  not  the  only  small-scale  feature  of  the  flow,  however.  As  noted 
by  Sandham  and  Yee  [20],  the  flow  also  develops  a  very  thin  region  of  high  strain  near  the 
stagnation  points. 

Although  there  is  a  physical  viscosity  in  the  flow  (in  this  case  the  Reynolds  number 
Re  =  100),  the  thicknesses  of  the  shock  and/or  high  strain  regions  may  eventually  become 
too  small  for  a  central  difference  scheme  to  handle.  Such  is  the  case  here  even  on  a  1024  X  192 
mesh.  In  the  presence  of  unresolved  gradients  central- difference  schemes  develop  oscillations 
which  lead  to  negative  temperatures  and  an  abrupt  halt  to  the  calculation. 

Indeed,  computations  with  both  the  spectral  and  the  sixth-order  schemes  on  32*,  64* 
and  128*  grids  develop  severe  oscillations  and  come  to  a  crashing  halt  between  t  —  105 
and  t  =  110.  Somewhat  curiously,  as  noted  by  Sandham  (1990,  private  communication), 
for  both  methods  the  oscillations  develop  first  not  in  the  vicinity  of  the  shock  but  in  the 
region  of  high  strain.  This  is  particularly  apparent  in  Figure  7,  which  shows  the  evolution  of 
the  pressure  for  the  spectral  method  calculation  (in  which  a  128*  grid  was  used  from  t  =  0 
to  i  =  75,  a  256  x  128  grid  from  t  =  75  to  <  =  100,  a  384  x  128  grid  from  t  =  100  to 
t  =  106.25,  a  768  x  144  grid  from  t  =  106.25  to  t  =  109.375,  and  a  1024  x  162  grid  from 
t  =  109.375  to  t  =  112.0).  The  computation  halted  shortly  after  i  =  112  due  to  negative 
temperatures  caused  by  severe  oscillations.  Figure  8  is  a  blow-up  of  the  central  region  at 
<  =  112.  Note  that  there  are  no  oscillations  apparent  in  the  vicinity  of  the  shocks.  Note  also 
that  the  oscillations  are  predominantly  in  the  x  direction.  Indeed,  an  examination  of  the 
spectral  coefficients  reveals  that  the  y  direction  is  quite  well-resolved.  The  price  of  resolving 
these  regions  with  a  non-dissipative  central  difference  scheme  can  easily  be  excessive,  as  the 
present  case  indicates. 
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The  vorticity  thickness  of  the  ENO  and  sixth-order  compact  methods  are  provided  in 
Figures  9  and  10,  respectively.  The  results  for  ENO  method  demonstrate  convergence  -  the 
vorticity  thickness  on  the  128^  grid  is  virtually  coincident  with  the  compact  method  result. 
The  errors  for  the  ENO  method  for  this  Mach  0.9  ccise  are  substantially  larger  than  for 
the  Mach  0.5  case  (see  Figure  1).  However,  this  is  due  primarily  to  the  diirerence  in  the 
stretching  parameter.  For  the  Mach  0.9  ENO  calculations,  a  weaker  transverse  stretching 
was  used  to  afford  finer  resolution  in  the  vicinity  of  the  shock. 

Results  for  the  lowest  4  Fourier  harmonics  for  the  two  methods  are  given  in  Figures  11 
and  12.  The  performance  of  the  two  methods  for  this  diagnostic  mimic  that  for  the  vorticity 
thickness. 

All  of  the  results  reported  thus  far  have  been  for  the  ENO  method  using  the  characteristic 
decomposition.  As  noted  at  the  end  of  Section  3,  componentwise  ENO  is  simpler  to  program 
and  is  less  expensive.  Figures  13  and  14  compare  the  two  versions  at  t  =  125  and  t  —  150, 
respectively.  The  componentwise  results  sutler  in  two  respects.  First,  the  shock  is  more 
diffuse.  In  fact,  at  t  =  125  the  shock  is  barely  visible.  Second,  there  are  appreciable  spurious 
oscillations.  Their  character  is  quite  different  from  the  oscillations  which  afflict  the  compact 
and  spectral  results.  They  are  far  less  regular  but  are  held  in  check  by  the  nonlinearly 
stable  adaptive  stencil.  Nevertheless,  the  small-scale  flow  features  of  the  componentwise 
ENO  results  are  quite  unreliable.  One  must  hesitate  to  use  this  method  for  applications 
in  which  the  small  scale  features  are  of  particular  interest,  such  as  transition  to  turbulence 
problems. 

A  comparison  of  the  characteristicwise  ENO  results  at  t  =  125  with  those  of  the  compact 
scheme  (Figure  6)  reveals  that  even  on  a  128^  grid  the  ENO  method  produces  a  numerical 
shock  thickness  which  is  much  larger  than  the  actual  thickness  for  this  viscous  problem. 
Moreover,  the  absence  of  a  shock  on  the  64^  grid  appears  due  to  the  delayed  flow  evolution 
(presumably  caused  by  the  inherent  viscosity  of  the  ENO  method)  that  is  apparent  in  Figure 
9. 


We  conclude  the  results  for  this  problem  with  Figure  15,  which  shows  the  long  time 
evolution  of  a  64^  characteristicwise  ENO  calculation  based  on  the  Euler  equations  (but 
starting  with  the  same  initial  conditions  as  the  Navior-Stokes  calculation  above).  Even  on 
this  coarse  grid,  and  without  the  aid  of  any  physical  viscosity,  the  ENO  method  exhibits 
solid  shock-capturing  behavior.  The  numerical  solution  shows  no  sign  of  nonlinear  instability 
and  spurious  small-scale  oscillations  are  absent. 


5  Isotropic  Turbulence 


In  [15],  Passot  and  Pouquet  simulated  two-dimensional,  compressible,  isotropic  flows  in  the 
turbulence  regime  using  a  Fourier  spectral  collocation  method.  They  identified  three  basic 
regimes:  shock-free,  weak  shocks,  and  strong  shocks.  Subsequently,  Erlebacher,  et  al.  [5] 
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developed  a  theory  of  compressible  turbulence  that  contained  a  more  refined  characterization 
of  the  different  regimes  of  compressible  turbulence  and  contained  a  useful  parametrization  of 
initial  conditions  that  permitted  precise  predictions  of  the  asymptotic  turbulence  state.  The 
direct  simulations  performed  in  these  studies  were  limited  to  quite  low  Reynolds  numbers, 
particularly  in  the  shock  regimes.  Gibbs  oscillations  arose  whenever  the  shocks  were  too  thin 
for  the  grid  to  resolve. 

Here  we  perform  simulations  of  compressible  isotropic  turbulence  using  both  spectral 
and  ENO  methods.  The  boundary  conditions  are  periodic  in  all  directions,  the  velocity 
is  normalized  by  its  initial  root-mean-square  value,  the  density  and  the  temperature  are 
normalized  by  their  mean  values,  the  pressure  as  for  the  free  shear  layer  problem,  and 
the  viscosity  is  taken  to  be  constant  /t  =  1/150.  We  compute  a  low  Mach  number  case 
where  the  shocks  are  weak  and  the  spectral  method  can  resolve  the  full  structure  with  256^ 
points.  Comparisons  between  different  ENO  schemes  and  between  ENO  schemes  and  spectral 
methods  are  furnished  for  both  large-scale  and  small-scale  flow  features 

In  Figure  16,  we  show  the  density  and  vorticity  contours  at  t  =  1  computed  with  the 
spectral  method  using  256^  points.  This  can  be  considered  to  be  a  resolved  solution.  Still, 
the  vorticity,  which  involves  derivatives  for  the  numerical  solution,  shows  some  oscillations. 
In  Figures  17  and  18,  we  show  the  density  and  vorticity  contours  at  t  =  1  for  the  spectral 
method  and  the  third  order  characteristic  ENO-LF,  respectively,  using  64^  and  128*  points. 
We  can  see  that  the  ENO  scheme  produces  non-oscillatory  results  but  the  spectral  method 
gives  noticeable  oscillations.  For  this  example,  the  component  ENO-com  produces  results 
similar  to  those  of  ENO-LF. 

In  Figures  19,  20  and  21,  we  show  the  time  history  of  the  average  Mach  number,  the 
maximum  Mach  number,  and  the  mininum  divergence,  for  the  spectral  schemes  and  the  third 
order  ENO-LF.  We  can  see  the  convergence  of  ENO  schmes  for  the  former  two  but  not  the 
latter. 

In  fact  minimum  divergence  is  achieved  exactly  inside  the  transition  regions:  to  resolve 
it,  one  has  to  resolve  the  full  transition  regions.  The  idea  of  using  high  order  shock  capturing 
methods  is  to  resolve  essential  features  in  the  smooth  part  of  the  flow  without  fully  resolving 
the  transition  regions  or  shocks.  An  important  topic  of  numerical  tests  is  to  verify  whether 
this  is  achieved.  For  this  example  we  indeed  achieve  this  as  indicated  by  Figures  19  and  20. 

We  then  compute  a  three-dimensional  version  of  this  problem.  We  present  in  Figure 
22,23,24  the  time  history  of  average  Mach  number,  maximum  Mach  number  and  minimum 
divergence  of  velocity.  These  have  characteristics  similar  to  their  two-dimensional  coun¬ 
terparts.  The  minimum  divergence  time  history  strongly  suggests  the  presence  of  three- 
dimensional  shocks.  However,  the  grid  resolution  is  not  sufficient  to  clearly  bring  them  out 
by  simple  flow  visualization. 


15 


6  Shock  Interaction  With  Entropy  Waves 


For  shock  interaction  with  weak  entropy  and  vorticity  waves,  some  qualitative  pictures  have 
already  been  presented  in  [23,  14],  which  are  similar  to  those  obtained  by  shock-fitting 
methods  [27].  Here  we  want  to  do  some  quantitative  studies  of  wave  amplification  factors. 
This  is  relevant  to  the  issue  raised  in  Example  1:  if  the  shock  or  the  rapid  transition  region  is 
not  completely  resolved,  can  we  still  resolve  smooth  information  passing  through  the  shock, 
such  as  the  amplification  factor  when  waves  pass  through  a  shock.  For  this  purpose  we  take 
an  entropy  wave  with  a  small  amplitude  so  that  the  linear  effects  dominate  and  a  comparison 
with  linear  amplification  factor  can  be  made.  For  third-order  ENO  with  150  x  20  grid  points 
(about  20  points  per  wavelength),  we  can  already  resolve  the  amplitude  amplification  factor 
to  within  5%.  This  compares  very  well  against  second  order  MUSCL  type  TVD  scheme  with 
the  same  number  of  grid  points  which  can  only  resolve  the  amplitude  amplification  factor 
with  an  error  six  to  ten  times  as  big.  Similar  comparisons  in  the  one-dimensional  case  were 
also  made  (via  graphs)  in  [23].  We  remark  that  this  quantitative  comparison  is  important  in 
this  case,  since  a  major  difference  between  ENO  and  TVD  schemes  is  that  the  latter  “clips” 
the  critical  points.  If  we  only  compare  contours  we  will  not  see  such  sharp  differences. 


The  details  of  this  problem  are  as  follows.  For  a  pure  shock  with  Mach  number  M  moving 
to  the  right,  we  add  an  entropy  wave 


P  = 


prC 


~^CO$Pr 


(34) 


where  =  kr{x  cos  -f-  y  sina,),  to  the  density  field  at  the  right  of  the  shock.  Here  a,  is 
the  angle  of  the  vortici'y  wave  with  the  shock,  kr  controls  the  number  of  waves,  and  e,  is 
the  scaled  amplitude.  In  order  to  enforce  periodic  boundary  conditions  in  the  y  direction, 
we  take  the  computational  domain  to  be  [0, 1]  x  [0, 


The  first  phase  of  this  computation  is  aimed  at  reducing  the  transients  that  arise  from 
the  discrete  ENO  approximation  to  the  moving  shock  wave.  We  run  the  scheme  until  the 
shock  moves  from  x  =  0.2  to  aj  =  0.8,  then  shift  the  data  leftwards  so  that  the  shock  is  again 
located  at  a  =  0.2,  and  repeat  this  process  six  times.  Then  for  each  fixed  x  to  the  left  of 
the  shock,  we  perform  a  Fourier  analysis  on  the  entropy  to  find  the  amplitude  e/,  where  (34) 
with  the  the  subscripts  “r”  replaced  by  "1”  denotes  the  entropy  wave  to  the  left  of  the  shock. 
The  resulting  amplitude  e/  is  then  averaged  over  an  x-interval  between  the  wave  front  and 
the  shock,  with  a  length  at  least  one  full  wavelength. 


The  computed  amplification  factors  f-,  together  with  the  linear  prediction  results,  for 
Mach  3,  Or  =  30°,  Cr  =  0.02,  kr  =  15,  are  listed  in  Table  4.  The  6th-order  ENO  method 
refers  to  a  scheme  which  is  sixth-order  in  space,  and  third-order  (Runge-Kutta)  in  time, 
with  a  reduced  CFL  number  of  0.2.  We  can  observe  from  Table  4  that,  especially  for  coarse 
grids,  higher-order  methods  indeed  produce  much  more  accurate  amplification  factors  than 
low-order  methods..  It  seems,  however,  that  we  cannot  reduce  the  error  below  a  certain 
threshhold  around  2%.  Again,  round-of  effects  might  be  playing  a  role  since  the  amplitude 
of  the  wave  is  far  smaller  than  the  shock  strength. 
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method 

N^  =  m 

=  100 

=  150 

2nd-order  MUSCL 

numH^i 

3rd-order  ENO 

6th-order  ENO 

3rd-order  ENO-2 

-47% 

■1.82% 

6th-order  ENO-2 

-8.43% 

-6.96% 

-2.38% 

Table  4:  Relative  Errors  in  Amplification  Factors 

7  Concluding  remarks 

ENO  schemes  based  on  fluxes  and  Runge-Kutta  type  TVD  time  discretizations,  introduced 
in  [22,  23]  are  implemented  on  Cray  2  supercomputers.  Vectorization  is  realized  for  all 
inner  loops.  Currently  the  code  runs  4.5  times  slower  than  the  classical  centered  difference 
schemes  of  the  same  order:  a  factor  of  2  is  due  to  the  upwind  flux  splitting  f  =  f'^  +  f~ , 
another  factor  of  2  is  due  to  the  adaptive  stencil  process.  If  characteristic  decompositions  are 
used,  the  CPU  time  is  increased  by  another  factor  of  1.7  to  3.  General  geometry  is  handled 
by  transformations.  Numerical  examples  include  2D  and  3D  homogeneous  turbulence,  shear 
flows,  and  shock  interaction  with  vorticity  waves.  ENO  schemes  show  their  advantage  when 
the  solution  contains  both  strong  shocks  and  detailed  structures:  with  a  relatively  coarse 
grid,  where  shocks  or  rapid  transition  regions  are  not  fully  resolved,  quantities  like  mini¬ 
mum  divergence  cannot  be  resolved,  but  the  numerical  result  is  still  stable  and  large-scale 
quantities  such  as  Mach  number  and  amplification  factors  can  be  well  resolved. 
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Figure  3.  Evolution  of  the  4  lowest  Fourier  harmonics  for  the  Mach  0.5  free  shear  layer 
problem  using  the  3rd-order  ENO  scheme. 
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Figure  4.  Evolution  of  the  4  lowest  Fourier  harmonics  for  the  Mach  0.5  free  shear  layer 
problem  using  the  6th-order  compact  scheme. 
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Figure  7.  Pressure  field  for  the  Mach  0.9  free  shear  layer  problem  using  the  spectral  scheme. 


Figure  8.  Blow-up  of  the  central  region  from  the  spectral  results  for  the  free  shear  layer 
problem  at  i  =  112. 
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Figure  11.  Evolution  of  the  4  lowest  Fourier  harmonics  for  the  Mach  0.9  free  shear  layer 
problem  using  the  3rd-order  ENO  scheme. 
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Figure  12.  Evolution  of  the  4  lowest  Fourier  heirmonics  for  the  Mach  0.9  free  shear  layer 
problem  using  the  6th-order  compact  scheme. 
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for  the  Mach  0.9  free  shear  layer  problem  at  t  =  125  using  the 
characteristic  form  (top)  and  in  component  form  (bottom).  Both 
results  are  given. 


64  x  64  128x128 
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left)  and  128^  (right)  results  are  given. 


t  =  100  ht  =  125 


Figure  15.  Pressure  field  for  the  Mach  0.9  free  shear  layer  problem  from  a  64*  inviscid  ENO-2 
calculation. 


Figures  16;  density  (left)  and  vorticity  (right)  contours  for  the  spectral  scheme  with  256^ 
grid  points. 
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Figures  17:  Density  (left)  and  vorticity  (right)  contours  for  the  spectral  scheme  with  64^ 
(top)  and  128*  (bottom)  points. 
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2-D  Compressible  T urbulence 


Figure  19:  Time  history  of  average  Mach  number. 
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Figure  20:  Time  history  of  maximum  Mach  number. 
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Figure  21:  Time  history  of  minimum  divergence. 
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Figure  24:  Time  history  of  minimum  divergence,  3D. 
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